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Abstract 

The reduced-density-matrix method is a promising candidate for the next generation electronic 
structure calculation method; it is equivalent to solve the Schrodinger equation for the ground 
state. The number of variables is the same as a four electron system and constant regardless of the 
electrons in the system. Thus many researchers have been dreaming of this much simpler method 
for quantum mechanics. In this chapter, we give a overview of the reduced-density-matrix method; 
details of the theories, methods, history, and some new computational results. Typically, the 
results are comparable to the CCSD(T) which is a sophisticated traditional approach in quantum 
chemistry. 
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I. INTRODUCTION 



Chemistry is an important branch of science which treats change of matter. It explains, 
for example, why and how a protein works, the process C0 2 converts to 2 , etc. The goal 
is to predict, understand, and control what will happen when we mix some substances. 
To do that, we usually do experiments which can be explosive, poisonous, expensive and 
unstable. Therefore, it is desirable to do chemistry without experiments. Fortunately, the 
basic equation for chemistry is already known and it is called the Schrodinger equation [12]. 
It is possible to approximately solve the Schrodinger equation using computers with various 
methods and algorithms. Such branch of chemistry is called quantum chemistry and it is 
our main interest [68]. 

Determining the exact or approximate solution to the Schrodinger equation is the funda- 
mental problem in quantum chemistry. This solution is called the wavefunction, or some- 
times referred as electronic structure. If we know the electronic structure, we can do chem- 
istry. Often such methods are referred as ab initio (Latin word which means "from the 
beginning") or the first principle method if approximations are not heuristic or do not em- 
ploy parameters from experiments. 

The ground state energy calculation of a non-relativistic and time-independent, TV- 
electron molecular system under the Born-Oppenheimer approximation is the most impor- 
tant problem [57]. It can be obtained as the lowest eigenvalue E of the electronic Schrodinger 
equation: 



tron % and nucleus A, and r^ is the distance between two distinct electrons. The solu- 
tion of (1), *(z) in L 2 (K N ), K = R 3 x \} with the inner product * 2 («)) = 
J ^i(z)^ 2 (z)dz, z = (x,s) e K (J dz includes integration over spin variables), is the 
wavefunction, and the corresponding eigenvalue E, the total energy of the system. 



HV(z) = EV(z), 



(1) 



where H is the Schrodinger operator or Hamiltonian defined by 





in which Za is the atomic number of the nucleus A, riA is the distance between the elec- 
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Besides, since electrons are fermions, the wavefunction itself is antisymmetric due to Pauli 
exclusion principle: 

. . ■ , Zi, . . . , Zj, . . . , Z N ) = —V(Zi, ...,Zj,...,Zi,..., z N ). 

That is to say, we must solve the Schrodinger equation in the antisymmetric subspace of 
L 2 {K N ). We denote such space as AL 2 (K N ). 

Even on computers, treating the iV-particle wavefunction is very difficult. Thus, we 
discretize the Hilbert space AL 2 (K N ) by taking antisymmetric products of the one-particle 
Hilbert space L 2 (K), whose complete orthonormal system (CONS) is {ipi\ c *l 1 . Each ipi is 
called single- electron wavefunctions or spin orbitals 

^:K->R. (i = 1,2, ...,oo) (3) 

We can explicitly construct a CONS of AL 2 (K N ) using {ipi}^! by the Slater determinants 
defined as follows [57]: 

1ph{Zl) A 2 ( Z l) ^jv( 2 l) 

V'iiM ^i 2 (z 2 ) ■■■ ^i N (z 2 ) 

Ipi^ZN) i) i2 {z N ) ■■■ i) iN {z N ) 

Here, we used an ordered set of indices / = {ix, . . . ,ij, . . . ,i k , . . . ,i N } where ij < i k , ij, ik € 
N. It is known that {^/} is a CONS of AL 2 (K N ) [33]. 

A second approximation to solve the Schrodinger equation would be choose carefully r 
functions from a CONS {ipi}^ by chemical or physical intuition. Then we construct a 
subspace of AL 2 (K N ) by the Slater determinants considering all possible combinations of N 
spin orbitals among r possibilities. In this case, solving the Schrodinger equation becomes the 
eigenvalue problem of the Hamiltonian matrix which now seems to be feasible on computers. 
Nevertheless, the dimension of the problem becomes r\/(N\(r — N)\) which obviously is 
impractical even for small values. The (approximate) ground state energy obtained by this 
procedure is considered the reference value, and called Full Configuration Iteration (Full CI) 
(energy). The mainstream approaches in quantum chemistry can be roughly interpreted as 
linear or nonlinear approximations of this eigenvalue problem, e.g., Hartree-Fock method, 
second-order perturbation methods, coupled cluster methods, truncated CI methods, etc. 
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Our main motivation is employ the second-order reduced density matrix (2-RDM) as the 
basic variable for quantum mechanics to construct simpler methods. Since only two-body 
interactions exist in nature, we can calculate all observables using the 2-RDM. Moreover, 
the number of variables of the 2-RDM is always four, regardless of the number of electrons 
in the system (r 4 when discretized), whereas the wavefunction scales like N (r\/(N\(r — N)\ 
when discretized). 

This chapter is organized as follows. In Section 2, we define the first-order and second- 
order reduced density matrices, and introduce the notion of iV-representability and its con- 
ditions. The reduced-density-matrix method is a viable implementation to approximate the 
ground state energy of molecular systems. The reduced-density-matrix method is formu- 
lated as an semi definite program in Section 3, and its numerical results using a parallel 
optimization code is given in Section 5. In Section 4, we give a brief historical note of this 
approach. Finally in Section 6, we give some concluding remarks. 

II. THE REDUCED-DENSITY-MATRIX METHOD 
A. Pure states and ensemble states 

Most generally, a quantum system containing N particles is described by the density 
matrix D which was introduced independently by von Neumann, Landau, and Bloch, and 
is an ensemble average of wavef unctions, 

D(z u z 2 ,..., z N , z[, z' 2 , ...,z' N ) = ^WiViiZu z 2 ,..., z N )^*(z[, z' 2 , z' N ), 

i 

where Wi > 0, Yli^i w i = 1> anc ^ is a CONS of an iV-particle state. For a pure state 

where the system is described by the wavefunction \1>, D can be written by 

D(zi, z 2 ,..., z N , z[, z' 2: . . . , z' N ) = ty(zi, z 2 ,..., z N )^*(z[, z' 2: . . . , z' N ). 

This is equivalent to requiring D to be idempotent; D 2 = D. Hereafter, when we refer to a 
state, it means be an ensemble if not otherwise specified. We are mainly interested in the 
pure state but usually we do not care about whether D is a pure or ensemble. This only 
becomes a problem when the system is degenerated or we consider its subsystems. 
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B. The first-order and second-order reduced density matrices 

Given an ensemble D(-), the first- order Reduced Density Matrix (1-RDM) [12] is defined 

by 

7O1, z[) = N J D(z u z 2 , . . . , z N , z[, z 2 ,..., z N )dz 2 dz 3 ■ ■ ■ dz N . 
The second-order Reduced Density Matrix (2-RDM) [28, 33, 35] is given by 

r(«i, z 2 , z ± , z 2 ) = / D(zi, z 2 , . . . , z N , z ± , z 2 , . . . , z N )dz 3 ■ ■ ■ dz N , 

and higher-order RDMs are defined in an analogous way. The normalization factor for 
the p-th order reduced density matrix is then p ^L p y ■ The second-quantized versions are 
defined using a set of creation and annihilation operators {a^aj}?^, where a\ creates and 
dj annihilates a particle at tpi(z) of D as follows [57]: 

7] = tr(a\a j D) = Y l w p (%\a\a j \%), 
v 

T Z = ^tr(a\a]a e a k D) = ^^w p (%\a\a]a e a k \V p ). 

v 

The normalization factor for the p-th order reduced density matrix is l/p\ for this case. The 
equivalence of these two different expressions can be found by 

7] = J $(z 1 )'y(zi,z' 1 )ij> j (z' 1 )dz 1 dz' 1 , 

T\\ = J ij*(z 1 )ij*( y z 2 )T(z 1 , z 2 , z[, z'^ t {z'^ k {z 2 )dz x dz 2 dz' x dz' 2 , 

where we used the single-particle wave function {-^(2)}°^ of (3). The following conditions 
are inherited by these definitions: 

(1) the 1-RDM and 2-RDM are Hermitian, 

7j = (7?r, rji = (rjfr, 

(2) the 2-RDM is antisymmetric, 

pu _ _pi* _ _p*i _ pi* 

k£ k$, $.k $.k ' 

(3) trace conditions are valid, 



(4) a partial trace condition holds between the 1-RDM and 2-RDM, 

k 

Additionally, we can find more conditions from the symmetry of the system. In particular, 
the spin symmetry of the (ground state) molecular systems is important and formulated as 
follows: 

(5) the total spin S 2 ; the 2-RDM should be the eigenstate of spin operator 

tv(S 2 T) = S(S + 1), 

where S 2 is defined as follows: 

S 2 = S 2 + S 2 + S 2 = S z + S 2 + S—S+ 




+ a l/3 a i* a ]a a jP 
ij 



The indices ia, j/3 means that we choose spin eigenfunctions of the z-axis for 
{fpiiz)}^, and reorder them so that % means i-th spacial function and a, /3 denote 
eigenfunctions of a-spin and /3-spin; {t/; ia (z), ^^(z)}^, respectively. 

(6) The z-component of the spin, S z can be chosen as integer or half integer, 

i 

In the subsequent discussion, one will notice that the 1-RDM can be disregarded throughout. 
However, we explicitly use it in order to prioritize the compactness of the notation. 

C. Solving the ground state problem using 1- and 2-RDMs 

The Hamiltonian of the most general form in second-quantization can be written by: 

H = ^2 v ) ai<x ) + 2 J2 w ke a i a A a l 

ij ijk£ 
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where v * and w^ e are one- and two-particle terms which can be calculated from the molecular 
Hamiltonian (2) by the Slater's rule as follows: 

where is the distance between an electron and a nucleus. Then, the total energy E can 
be expressed using D as follows: 

E = tr(HD). 

The ground state energy E min can be calculated minimizing the total energy over 1- and 
2-RDMs. 

E min = mintr(HD) 

= min tr(J2 v ) a A + \ Yl w ki a i a A a l) D ( 4 ) 

ij ijki 

= min^v]tr(a i a] J D) + Yw l ^tr(a i a j a\a\D) 

ij ijki 

= n>i"{£^i-E ? <«}- (5) 

It is easy to show that this minimization (4) is equivalent to solve the Schrodinger equation 
for the ground state, and often such kind of methods are called variational methods. More- 
over, (5) is also equivalent to solve the Schrodinger equation for the ground state. Advantage 
of using 1- and 2-RDMs instead of D is that the number of variables are reduced drastically. 

D. The iV-representability problem and the iV-representability conditions 

In 1950's, researchers based on the above facts chose the 1- and 2-RDMs as basic variables, 
and did some variational calculations; according to Lowdin [34], F. London, J. E. Mayer, 
A. J. Coleman, P. O. Lowdin, R. McWeeny, N. A. March, C. A. Coulson and others have 
attempt to minimize via (5). However, their results were considerably lower than the true 
energy. The reason is that the trial 2-RDMs were not actually derived from an existing 
density matrix D. We need some more conditions on the trial 2- RDM to ensure that it 
comes from a true D. Such formalism of the problem was first described by A. J. Coleman 
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in 1963 and named N -representability problem, and these conditions are known to be the 
N -representability conditions [9]; Given a trial p-th order RDM T^ p \ if there exists some 
wavefunction or ensemble which reduces to the p-th order RDM T^' p \ then this r( p ) is pure 
or ensemble iV-representable, respectively. 

E. On the complete iV-representability conditions 

Given a 1- and 2-RDM, they should satisfy the relations (1) to (4) of Section 2.2 The 
1- and 2-RDM for the ground state should additionally satisfy (5) and (6). Therefore, 
these conditions are necessary conditions for the iV-representability. Unfortunately, these 
conditions are not sufficient, thus early attempts failed and obtained very low energy. The 
necessary and sufficient condition for the 1-RDM is relatively easy [9, 30]. However, the 
complete (sufficient) iV-representability of the 2-RDM is very complicated in general. Garrod 
and Percus [26] showed that the 2-RDM T is ensemble iV-representable if and only if 

ijkl 

where H u is every possible Hamiltonian and -E^in ls the ground state energy corresponds to 
H". Thus, the ensemble iV-representable set £ N can be defined by: 

S N = {T | tr(ifT) > E^ in ,for all possible H" and E^J. 

This result is theoretical and very important, but totally not practical since if one wants to 
calculate the exact ground state energy of a Hamiltonian, then he/she must know the exact 
ground state energy of the system beforehand. This is a tautology! After that, many re- 
searchers seek the complete iV-representability condition, and did not succeed. A meaningful 
result from complexity theory was obtained by Liu et al. [32] in 2007. They showed that 
the computational complexity of the iV-representability problem is QMA-complete, which 
is the quantum generalization of NP-completeness. Thus it is almost hopeless to find an 
efficient algorithm to decide whether a given 2-RDM is iV-representable or not. We can 
consider a more physical example: the ground state problem of the spin-glass Hamiltonian 
is known to be a very hard problem, and equivalent to solve the max-cut problem or the 
traveling sales person problem, which in turn are known to be NP-hard [3]. If the complete 
7V-representability conditions were easy to handle, we could solve such difficult problems in 
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tr(H 1 r)>E 1 



A/-representable 
set 




trfH^E 2 



tr(H 3 r)>E 3 tr(HY)>E 4 

FIG. 1: Schematic representation of the iV-representable set; a Hamiltonian and its ground state 
energy serves as a characterization of the set. 

computer science as well. Currently we do not know how to solve these problems efficiently. 
We just want to stress, the complete 7V-representability is a really hard problem. 

F. Formulating as a variational problem, and its geometrical representation 

The problem we want to solve can be formulated using the 1- and 2-RDMs as basic 
variables, 

-Etain = min tr(HD) 

ij ijki 

£ N is known to be a compact convex set. Besides, all possible Hamiltonians and the 
corresponding ground state energies serves as a characterization of this convex set. To be 
more precise, a 2-RDM corresponding to the ground state of an iV-particle Hamiltonian 
is a surface point, and any surface point of £ N corresponds to the ground state of some 
Hamiltonian [53]. The compact and convex set of the iV-representable set is represented 
as an ellipse (we do not show, but there are also cusp points as well) in Fig. 1. The 
Hamiltonians H l , H 2 , H 3 , and H 4 , and their ground state energies E 1 , E 2 , E 3 , and E 4 
serves as iV-representability conditions, respectively. 
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G. Some known necessary iV-representability conditions 

We should not be demotivated by the facts of the previous subsection. Understanding the 
chemical and/or physical meaning of the necessary iV-representability conditions is much 
more important. Mathematical theorems do not tell about chemistry or physics. In practice, 
iV-representability conditions for the molecular systems might not be so difficult. 

We seek for chemically and/or physically meaningful necessary iV-representability con- 
ditions on the 1- and 2-RDMs. The necessary and sufficient conditions for an ensemble 
iV-representability for the 1-RDM is characterized by its eigenvalues lying between and 
1 [9, 30]. For the pure state, it is more complicated [1]. Coleman introduced two neces- 
sary conditions called the P and Q conditions [9]. These conditions require the positive 
semidefiniteness of the P-matrix (r), and the Q-matrix 

Pl{ = tr(a]a]a e a k D) y 0, (6) 

Qkt = H a i a A a l D ) °- ( 7 ) 

Another important necessary condition is called the G condition [26], which also require 
positive semidefiniteness of the G-matrix defined as follows 

G\\ = ti(a\ aj a\a k D) y 0. (8) 

In the original paper by Garrod-Percus, the definition of the G-matrix is non-linear: 

G\\ = tr((a\ aj - 7 ;)(ak - Y k )D) y 0, (9) 

but for a fixed particle state, these G-matrices share the same eigenvalues, since 7* can be 
replaced with Yjjf 12i a \ a i [16]- In Zhao et al [67], we can find explicitly formula of Tl and 
T2 conditions from Erdahl's survey paper [14]: 

( T1 )tem = tr((ala]ala n a m a e + a n a m a e a\a]al)D) y 0, (10) 
( T2 )foL = ^{{a\a]a k a ] n a m ai + ala m aia\ala k )D) y 0, (11) 

which are stronger conditions. An important property of these matrices is that Q, G, Tl 
and T2-matrices can be written only by linear combinations of the 2-RDM elements like 
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following: 

Qm = m - - (« + *H) + m + «) - zrz, 

g % L = sH - zrg, 

(T2)f m = A[j, k]A[m, n](±M 7 i + ~ 2 

where A is the antisymmetrizer operator acting on an arbitrary function f(i,j, k), 

A[i,j, k]f(i,j, k) = f(i,j, k) - f(i,k,j) - f(j,i,k) + f(j,k,i) + f(k,i,j) - f(k,j,i). 

For Tl and T2's cases, the 3-RDM terms cancel out. The T2' condition replaces the T2 
condition and is slightly strengthened by the addition of the one-particle operator [5, 40]. 

Other positive semidefmite type representability conditions are known such as the B and 
C. However, they are implied by the G condition [31]. 

We can extend these conditions to positive semidefmiteness of higher order RDMs. These 
extensions seems to be known for a long time. Erdahl and Jin [17] formulated the p-th order 
approximation to the iV-particle density matrix in terms of the semidefmiteness conditions 
on the p-th order RDMs, which is an generalization of the P, Q, G, Tl and T2' conditions. 

H. The reduced-density-matrix method 

We call as the reduced- density-matrix method, the variational method having the 2-RDM 
(and the 1-RDM) as the basic variable(s) restricted to some approximation £ N of the N- 
representability set £ N . It can be formulated as follows: 

E min = ...in {]>>;-; + 5>£r&}- (12) 

i j ijkl 

Among the possibilities, we usually consider the set obtained by imposing some necessary 
conditions for the iV-representability. The set £ N should satisfy the following properties. 

• satisfies some necessary conditions of ensemble iV-representability. 

• compact set, so that a linear functional (the Hamiltonian) has a minimum. 

• convex set, so that the solution would not be stuck into local minima. 
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• stringent, so that resultant 2-RDM should be physically or chemically meaningful. 

• computationally feasible and/or efficient. 

• completely general: since the form of the Hamiltonian is totally general, it is not only 
applicable to chemistry but also to physics. 

• ab initio: no empirical parameters. Currently very successful methods based on the 
density functional theories employ a lot of empirical parameters. 

We can find new iV-representability conditions from chemical or physical requirements 
satisfying the above properties by constructing a Hamiltonian and obtaining an upper bound 
to the ground state energy. Then we can add this as a new condition. These "cuts" may 
strengthen S N . 

Trivial iV-representability conditions with the P, Q, G, Tl and T2' conditions and every 
possible combination of the P, Q, G, Tl, T2 and T2' conditions, satisfy the above criteria. 
These variational energies have the following property: 

• If we add more necessary conditions, the calculated energy usually becomes better and 
would never become worse: 

EpQ < Epqg < EpQGTl < EpQGTlT2 < EpQQT\T2' < PfullCI, 

where Epq is the variational energy with the P and Q conditions, Epqgtit2 is the vari- 
ational energy with the P, Q, G, Tl and T2 conditions, etc. (see Fig 2). This property 
is totally opposite from traditional wavefunction approaches. Variational calculation using 
the wavefunction gives upper bounds, and approximation to the total energy becomes lower 
(better) when the variational space becomes larger . 

Note that the obtained 2-RDM may not be unique, even when the original problem is 
non-degenerated and the energy is unique. 

I. Some interpretations on conditions 

We usually enforce only the necessary conditions on the trial 2-RDMs. Therefore, the 
RDM method gives lower bounds to the exact energy, an iV-representable 1-RDM and a 



13 




tr(H 3 r)>E 3 tr(HV)>E 4 

FIG. 2: Schematic representation of approximate iV-representable sets; a better set shrinks. 

non-physical 2-RDMs. Thus it is important to realize the physical meaning of the necessary 
conditions employed in the calculations. The following results may be useful to interpret 
results from actual calculations on molecules and atoms: 

(a) If a trial 2-RDM V satisfies the P and Q conditions, the original 2-RDM is ensemble 

iV-representable [9]. 

(b) If a trial 2-RDM V satisfies the G condition, then 1-RDM from the original 2-RDM is 

ensemble iV-representable [42]. 

(c) If the Hamiltonian of a system is time-reversal invariant, and the number of particles 

N, is even, the necessary and sufficient condition that an approximate 1-RDM corre- 
sponding to a non-degenerate energy eigenstate be iV-representable is that its natural 
spin-orbital occupation numbers be equal in pairs [56] . Moreover Coleman proved that 
the AGP (anti-symmetrized power) wavefunctions covers all such 1-RDMs [10]. Thus 
if these conditions apply to the systems, we always obtain pure representable 1-RDMs 
with necessary iV-representability conditions. 

(d) The G condition is related to the AGP type wavefunction, and gives the correct energy 

for the Hamiltonian for which the ground state can be written by the AGP function 
[18]. Moreover, the AGP type wavefunction is closely related to the superconductivity 
[44]. 
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(e) The G condition is exact at the high-correlation limit for the Hubbard model, since the 
two particle term of the model (^ Xlf a !t a ^ a 4 a ^' where U > 0, L is the number of 
sites, i is the i-th site, and t and I denote up-spin and down-spin of the electrons, 
respectively) is a G-type Hamiltonian (Ylijke^ke 

), which is bounded by zero 

[64]. 

III. FORMULATION OF THE RDM METHOD'S PROBLEM AS AN SEMIDEF- 
INITE PROGRAM AND ITS SOLUTION BY THE INTERIOR-POINT METHOD 

A. Semidefinite program 

Semidefinite program (SDP) has established as an important class of problems in opti- 
mization since 1990's. It is known to have an elegant mathematical theory, and an efficient 
algorithm called interior-point method which can solve it in polynomial-time complexity. 
Refer for instance to [58] for a nice survey about SDPs. 

Let C, A p (p = l,2,...,w) be given block-diagonal real symmetric matrices with pre- 
scribed block sizes, b e M" 1 and c, a p e K s (p — 1, 2, . . . , s) be given real vectors. We denote 
by Diag(a) a diagonal matrix with the elements of the vector a on its diagonal. 

An SDP is defined for instance by 

f maximize tr(CX) + tr(Diag(c)Diag(cc)) 
subject to tr(A p X) + tr(Diag(a p )Diag(a3)) — b p , (p — 1, 2, . . . , u) (13) 
x y O, x E R s , 

where we refer it as the primal SDP. The notation X y O means that X is symmetric 
positive semidefinite. Then, we can define the dual SDP as 

( 

minimize b T y 

u 

subject to S = ^ A p y p — C y O, 

u P=1 (14) 
^Diag(a p )?/ P = Diag(c), 

P =i 

y G R u . 

The variables for the primal SDP is (X,x) while for the dual SDP is (S,y). Under mild 
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assumptions [69] [58], the solution of (13-14) should satisfy 

tr(A p X) + tr(Diag(a p )Diag(a;)) = b p , (p = 1, 2, . . . , u) 

u 

S = Y, A PVP-Ch0, 

p=l 

x,s y o, 

b T y - tr(CX) - tr(Diag(c)Diag(a;)) = 0. (15) 

These conditions are equivalent to earlier result by Erdahl [15] and Bellman and Fan [4]. 

The advantage of considering the variables in the primal and dual SDPs simultaneously 
is that we can check the numerical correctness of the approximate solution from the above 
relations. 

B. Formulation of the RDM method's problem as an SDP 

Hereafter, we assume that we have chosen r spin orbitals from (3), which is assumed to 
give a good approximation for the desired wavefunction we seek. There are plenty of these 
bases in quantum chemistry, and we actually use them on the numerical experiments which 
follows. Also, notice that all definitions and notions of iV-representability and its conditions 
can be defined accordingly using this finite basis of CONS. 

The RDM method's problem imposing some necessary iV-representability conditions such 
as the P (6), Q (7), G (8), Tl (10), T2 (11) or T2' conditions, is in fact an SDP. In order 
to make its formulation more clear, we perform some linear transformations on the matrices 
involved in the problem. In (12), the 1-RDM variable 7, and the corresponding Hamiltonian 
v have two indexes, which correspond to ordinary matrices in linear algebra. However, the 
other matrices involved in the calculations have four or even six indices each. To convert 
from these notations convenient for quantum chemists to the notation of elementary linear 
algebra, we need to map each pair or triple k) of indices to a composite index on 
these matrices. For instance, the 2- RDM element (1 < % < j ' < r; 1 < k < £ < r) will 
be mapped to ^ j-i+(2r-i)(i-i)/2,e-k+(2r-k)(k-i)/2, which results in a symmetric matrix of size 
r(r — l)/2 x r(r — l)/2. We assume henceforth that all matrices are transformed to become 
two-index matrices, and we keep the same notation as before for simplicity. Furthermore, 
due to spin symmetry [67], all these matrices will reduce to block-diagonal matrices of sizes 
specified in Table I [22, 23] [70]. 
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TABLE I: Sizes of the SDP as a function of the number of spin orbitals r for each necessary 
TV-representability condition. 

iV-reprcs. cond. size of block matrices 

7^0 r/2 X r/2 (2 blocks) 

I y 7 r/2 X r/2 (2 blocks) 

P condition (r/2) 2 X (r/2) 2 (1 block), f ^ J X I l- 1,1 " ks; 

Q condition (r/2) 2 X (r/2) 2 (1 block), j X | |2 :■: ■■ : • 

G condition 2 (r/2) 2 X 2 (r/2) 2 (1 block), (r/2) 2 X (r/2) 2 (2 blocks) 

Tl condition § f r ^ 2 J X § j (2 blocks), ^ 2 j X ^ 2 J (2 blocks) 

/ 3r/2 \ / 3r/2 \ / r/2 \ / r/2 \ 

T2 condition g I I X g I ' I (2 blocks), § I J X § I I (2 blocks) 

Tl' condition j + |j ^ 3r ^ 2 ^ x § + § ^ 3r ^ 2 j ( 2 block s), § ^ 2 j X ^ ^ 2 ] 2 i- ■ 

, . 2 /4 + 1 \ / r (r/2 - 1) /4 + 1 \ 



s in 



(14) 5 + 2 



r/2 + 1 
2 



here | j — j^7^jTb\ > f° r integers a > 6 > 0. 



Now, let us define a linear transformation svec: 

g„ _^ R v{v+i)/2 from the 

space of v x t> 

symmetric matrices S 1 ". For U EE> V , 

svec(C7) = (C/ n , a/2«/i 2 , ^22, ^33, • • • , V2U lv , U VV ) T . 

Then, defining y = (svec(7) T , svev(r) T ) T e R u , b = (svec(v) T ,svev(w) T ) T G R", and 
finding the suitable matrices C, A p {p — 1,2, ... ,u) and vectors c, a p {p = 1,2, ... ,u) for 
the corresponding necessary iV-representability conditions of Table I, for instance, we can 
cast the problem (12) as (14). 

Although these transformations and formulations seem a little confusing, it is in fact the 
formulation Garrod et al. arrived 35 years ago [24]. Nakata et al. [46-48] formulated the 
problem as the primal SDP (13) instead when they performed the first computation as an 
SDP. For a more detailed discussion about these transformations and the formulations, see 
for instance [22, 23, 67]. 



C. Theoretical computational complexity of the primal-dual interior-point method 

As it was previously mentioned, SDPs can be solved in polynomial-time using interior- 
point methods [58]. In particular, employing the parallel code SDPARA [20], which is 
an implementation of the primal-dual interior-point method, one can theoretically expect 
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that it will take 0( y /v milx \oge~ 1 ) iteration with 0(u 2 f 2 /d + u 3 /d + + t> r 3 nax ) floating- 

point operations per iteration. Here t> max refers to the size of the largest block matrix in 
A p (p — 1, 2, . . . , u), f is the maximum number of nonzero elements in each of these matrices, 
d is the total number of available CPU cores in the parallel computer, and e is the accuracy 
which we can expect when we replace the rhs of (15) "=0" by "< e" (where e > 0). In our 
case, u = 0(r A ), t> max = 0(r 3 ), f = 0(1), and therefore, the total theoretical floating-point 
operations is £>(r 13 - 5 loge^/d) [23]. 

IV. SOME HISTORICAL REMARKS 

Here we make an attempt to list some articles related to our work in chronological order. 
However, it is far from being complete. 

The definition of the RDM was explicitly spelled out by Husimi [28] in 1940. The depen- 
dence of the energy on the 2-RDM (and 1-RDM) appeared in Lowdin [33] and Mayer [35] 
in 1955. The necessary and sufficient conditions for an ensemble 1-RDM, i.e., ^ 7 ^ I, 
tr(7) = N, were obtained by Kuhn [30] in 1960 and Coleman [9] in 1963. In this latter 
article, the precise formulation of the iV-representability problem, the P and Q conditions 
(for the 2-RDM) were also stated. In the next year, the G condition was proposed by Garrod 
and Percus [26]. 

The restriction of the iV-representability problem only on the diagonal elements of the 
2-RDM, known as the diagonal problem, were investigated by Weinhold and Wilson [63], 
Davidson [11], McRae-Davidson [41], and Yoseloff [66] since late 1960's. Other progresses 
on this topic can be found in the survey [2]. 

The first variational calculation on the 2-RDM imposing the necessary iV-representability 
conditions for the doubly ionized carbon C ++ (N = 4) were performed by Kijewski and co- 
authors since late 1960 's [29] (see earlier reference therein). Garrod and co-authors proposed 
several algorithms, some of them which resemble modern optimization algorithms, and re- 
ported results for the beryllium (N = 4) [24, 25, 54]. In particular, Mihailovic and Rosina 
applied it to nuclear physics [43] , but obtained large deviations to the full CI calculations if 
compared to electronic systems. 

In the 1978 survey paper of Erdahl [14], we can find the conditions which is knows as 
Tl, T2 [67], and T2' conditions [5, 40]. He also proposed algorithms based on the exact 
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mathematical characterization of solutions of the lower bound method (RDM method) in 
the next year [15]. 

These were the golden ages for the RDM research, but somehow faded away because it 
was realized soon that the underlying problem is computationally difficult and poor results 
were obtained for nuclear systems. 

A revival of the 2-RDM approach appeared since 1992 due to Valdemoro [59], Nakatsuji 
and Yasuda [51], and Mazziotti [36]. This approach is based on the density equation or 
the contracted Schrodinger equation (CSE), which is equivalent to solve the Schrodinger 
equation. Nakatsuji proved that if an iV-representable 4- RDM satisfies the CSE, then the 
original D satisfies the Schrodinger equation and vice versa [50]. Valdemoro, Nakatsuji- 
Yasuda and Mazziotti consider the 2-RDM as the basic variable. The CSE requires 1- to 
4-RDMs, thus they reconstruct 3- and 4-RDMs using 1- and 2-RDMs and solve the CSE 
iteratively. In this approach, they assume the resultant 2-RDM is nearly iV-representable 
because the reconstruction functional is physically relevant; they did not explicitly impose 
any iV-representability conditions. Their results are quite good, and can be compared to 
single and double CI for small atoms and molecules such as Be, Ne, and CH 3 F. The absolute 
values of negative eigenvalues of the P, Q, and G-matrices were small. However, researchers 
payed little attention to this approach since non-convergence or divergence occur especially 
where the correlations are strong [13]. There are difficulties in systematic refinements of the 
reconstruction functional even though some improvements are reported [65], but there are 
even more miscellaneous problems [52]. 

In 2001, Nakata et at. [48] were the first to employ an optimization software to solve the 
RDM method's problem as an SDP, and reported computational results imposing the P, Q, 
and G conditions on a series of small atoms and molecules. The results were better than 
the SDCI calculations, and obtained correlation energies from 100% to 120%. As mentioned 
in Section III A, these results have a numerical certificate of correctness, which could not 
be obtained before the advent of interior-point methods. Also, the numerical convergence 
does not depend on the initial guess as it is common in Hartree-Fock, CCSD methods, etc. 
Moreover, there exists a global minimum, whereas this is not guaranteed in the CSE approach 
for instance. In the next year, Nakata et al. [46] showed results for the dissociation limit 
for several molecules including triple bonded N2, demonstrating numerically that the RDM 
method do not break down as the single reference methods such as CCSD and perturbation 
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methods. However, it was also shown that size consistency is slightly deviated. 

The inclusion of the Weinhold- Wilson inequalities [63], which was not satisfied only in- 
cluding the P, Q, and G conditions [46], however showed little progress in the results [47]. 

In 2002, Mazziotti immediately reproduced Nakata et a/.'s results and applied to diatomic 
molecules [37]. The prolific research by Mazziotti and his colleagues in the following years 
[40] corroborated with these results. 

In 2004, a breakthrough was done by Zhao et al. They included additionally the Tl 
and T2 conditions, which became very strong conditions for small molecules and atoms. 
They noted a "spectacular increase in accuracy" and results were comparable to CCSD(T); 
typically the correlation energies for various atoms and molecules were between 100% to 
101% [22, 45, 67]. 

In the same year, Mazziotti announced the RRSDP method [38] in which he reformu- 
lates the SDP problem as an nonlinear and nonconvex problem and applies a quasi-Newton 
method to solve it [6]. In 2006, Cances et al. proposed and implemented the dual problem 
of (12) [7]. 

Applications to the one-dimensional Hubbard model was done by Hammond et al. [27]. 
They calculated the Hubbard models with P, Q, G, and T2 conditions up to 14 sites. The 
obtained error per site was —0.0089 for L = 14's case with P, Q, G, and T2 conditions for 
U — 8, when correlation is strongest. Nakata et al. investigated the high correlation limit 
using multiple-precision arithmetic version of SDP solver, called the SDPA-GMP [45]. At 
the high correlation limit, they reproduced the exact energy and proved that the G condition 
is exact [64] . 

The size-consistency and size-extensivity are important properties when the size of the 
systems becomes bigger or larger. Nakata et al. found slight deviations [46], but Van 
Aggelen et al. showed very clear and systematic examples that the RDM method gives 
incorrect dissociation limit with fractional charges on the well-separated atoms of diatomic 
molecules with P, Q and G conditions. For the CO's case, the Mulliken populations were 
5.98 and 8.00 at the dissociation limit. Even adding Tl and T2 conditions they did not fix 
the problem [60, 61]. Nakata and Yasuda investigated numerically that size-extensivity is 
also slightly violated by calculating 32 non-interacting CH 4 and N 2 [49]. The inextensive 
contributions to energies are 3 x 10~ 4 and 3 x 10~ 3 a.u. using the STO-6G basis set, 
respectively. Later, Verstichel et al. also proposed a method to "cure" this pathological 
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behaviour of the RDM method, however quite demanding [62]. Currently solutions to size- 
consistency is not practical. 

V. NUMERICAL RESULTS FOR THE RDM METHOD 

Here, we give numerical results we obtained so far for the RDM method imposing some 
necessary iV-representability conditions. Some of them are completely new. 

A. New numerical results for larger systems 

We present here some numerical results for the largest systems solved so far by our group. 

The SDPs obtained by the RDM method imposing the P, Q, G or P, Q, G, Tl, T2' 
conditions were solved using the parallel code SDPARA 7.3.2 [20, 21]. The calculations were 
performed at the Kyoto University's T2K supercomputer using 128 nodes, were each node 
has 4 CPUs (quad-core AMD Opteron 8356 2.3GHz) and 32GB of memory, giving a total of 
2048 cores; and at a self-made computer cluster using 16 nodes, were each node has 2 CPUs 
(quad-core Intel Xeon 5460 3.16GHz) and 48GB of memory, giving a total of 128 cores. 

Table II shows the results for five molecules were r = 28, 30 or 36 spin orbitals were 
used [21]. The full CI and SDCI (singly and doubly substituted configuration interaction) 
calculations were performed using the package Gamess [55], while CCSD(T) (coupled cluster 
singles and doubles with perturbational treatment of triples) and Hartree-Fock calculations 
were obtained by Gaussian98 [19]. The entries, excepting the full CI, give the ground state 
energy differences to the full CI. The RDM method always gives a energy lower than full 
CI, while SDCI and Hartree-Fock give higher. CCSD(T) usually results in higher energy, 
but not necessarily. Units are in Hartree. The acceptable accuracy in quantum chemistry is 
lkcal/mol which corresponds to approximately 0.0016 Hartree. Also, the correlation energy 
5 corr is an important measure in quantum chemistry. It is defined as a percentage relative 
to the Hartree-Fock (0%) and full CI (100%): 

E — -Ehf 

£corr = ~=, = X 100, 

— -&FCI 

where E the energy calculated by the RDM method, CCSD(T) or SDCI. 

From Table II, we can conclude that the RDM method imposing the P,Q,G,T1,T2' 
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TABLE II: Ground state energies (in differences from that of full CI) calculated by the RDM 
method imposing the P, Q, G, Tl, T2' conditions from SDPARA 7.3.2, and those obtained by 
CCSD(T), SDCI, and Hartree-Fock from Gamess and Gaussian98. The last column shows the full 
CI energies. The energy units are in Hartree (= 4.3598 x 10~ 18 J). The correlation energies (0% for 
Hartree-Fock and 100% for full CI) in percentage are also shown in the second row. 



system state 


basis 


r N(N a )2S + 


1 A£ ' PQGT1T2 1 


AE CCSD(T) 


AEsdCI AE HF 


Bfci 




NH 2 


'A, 


double-C 


28 10 (5) 


1 


-0.000 6 


+0.000 63 


+0.008 74+0.141 98 


-55.624 71 












100.4 


99.55 


93.84 


100 




CH 2 


'i, 


double-^ 


28 8 (4) 


1 


-0.000 4 
100.4 


+0.000 59 
99.42 


+0.005 80+0.100 67 
94.24 


-38.962 
100 


24 


NH 3 


l A, 


valence doublc-C30 10 (5) 


1 


-0.000 5 


+0.000 49 


+0.007 46+0.128 75 


-56.304 


89 












100.4 


99.62 


94.45 


100 




CH 3 


2 A 


valence double-C30 9 (5) 


2 


-0.000 3 


+0.000 31 


+0.004 01 +0.094 54 


-39.644 


14 












100.3 


99.67 


95.75 


100 




c 2 




valence double-^ 36 12 (6) 


1 


-0.003 5 


+0.000 39 


+0.055 98+0.285 66 


-75.642 


11 












101.2 


99.86 


80.41 


100 





conditions, gives equally better energies than CCSD(T) in absolute value. C 2 molecule is 
exceptional, however, it is known to be a difficult system in quantum chemistry. 

Table III shows the same result for the O^" molecule were r = 40 spin orbitals were used. 
The full CI calculation was not possible for this case due to the size limit on the computer, 
and therefore, we restrict to show only the ground state energy corresponding for each entry. 

TABLE III: Ground state energies calculated by the RDM method imposing the P, Q, G conditions 
from SDPARA 7.3.2, and those obtained by CCSD(T), SDCI, and Hartree-Fock from Gamess and 
Gaussian98. The energy units are in Hartree (= 4.3598 x 10 _18 J). 



system state basis r N(N a )2 


:S + 1 Ep QG 


E CCSD(T) E S DCI 


Ehf 


2 n g doublc-C40 15 (8) 


2 -149.450 2 


-149.385 95 -149.360 2 


6-149.091 83 



Table IV shows the typical sizes of the problem when formulated as an SDP (see Sec- 
tion III A), and the computational time to solve them by SDPARA 7.3.2 at the T2K super- 
computer or at the computer cluster. 

In the previous work [45], we only could solve the RDM method's problem with 
P,Q,G,T1,T2' conditions up to r = 28 spin orbitals. Here, we give results up to r = 36. 
This achievement was possible due to a major update in the parallel code SDPARA 7.3.2 
[20, 21]. It became faster, and now it can take advantage of multi-core (multi-thread) com- 
putation in addition to the ordinary MPI (message passing interface) computation. 

For other physical properties such as the dipole moments, refer to [45] . 
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TABLE IV: Size of SDPs obtained by the RDM method imposing some A-representability condi- 
tions shown at Tables II and III, and their computational time when solved at the T2K supercom- 
puter or at the computer cluster (c.c). 



system 


r A-repres. cond. 


U "max 


time (s) system CPU 


cores 


NH 2 " 


28 P,Q,G,T1,T2! 


27,888 4,032 


27,949 


T2K 


2048 


CH 2 


28 P,Q,G,T1,T2' 


27,888 4,032 


26,656 


T2K 


2048 


NH 3 


30 P,Q,G,T1,T2' 


36,795 4,965 


72,026 


T2K 


2048 


CH 3 


30 P,Q,G,Tl,T2> 


36,795 4,965 


68,593 


T2K 


2048 


c 2 


36 P,Q,G,T1,T2' 


76,554 8,604 1,554,675 


c.c. 


128 


0+ 


40 P, Q, G 


116,910 800 


5,943 


T2K 


2048 



B. Summary of the numerical experiments 

We present a graphical summary of the data obtained in our previous work [45] with the 
addition of new ones presented in the previous section. 

The ground state energy differences to full CI of the RDM method imposing the P,Q,G 
or P,Q,G,T1, and of Hartree-Fock for the 57 atomic or molecular systems [45] and those 
5 shown at Table II are plotted in Figure 3. Each horizontal bar corresponds to a system 
and they are ordered accordingly to the order it appears in the tables. That is, the lowest 
one corresponds to the Lithium atom with r = 10 spin orbitals [45], while the upper 5 
corresponds to the data of Table II (notice that we do not have the values for P, Q, G and 
P, Q, G, Tl entries for this case). 

Apparently, it does not seems to exist an correlation between the Hartree-Fock and the 
RDM method's results. However, we clearly notice that the imposing the P,Q,G conditions, 
we can obtain results much better than the Hartree-Fock's ones. 

Figure 4 shows the ground state energy differences to full CI of the RDM method imposing 
the P,Q,G,T1 (same value as Figure 3) or P,Q,G,T1,T2', and of CCSD(T) for the 57 
atomic or molecular systems [45] and those 5 shown at Table II (notice that we do not have 
the values for P, Q, G, Tl entries for this case). 

If you compare the RDM method with P,Q,G,T1,T2' conditions and the CCSD(T) 
values, they seems equally good. However, in some case CCSD(T) can fails to converge. 
There are 4 case in [45] which are replaced by zero in Figure 4. The largest deviation of 
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-0.05 0.05 0.1 0.15 0.2 0.25 0.3 

ground state energy differences in Hartree 



FIG. 3: Ground state energies (in differences from that of full CI) calculated by the RDM method 
imposing the P,Q,G, or P,Q,G,T1 conditions, and those obtained by Hartree-Fock for the 57 
atomic and molecular systems of [45] and the 5 of Table II. 




-0.03 -0.02 -0.01 0.01 0.02 0.03 

ground state energy differences in Hartree 

FIG. 4: Ground state energies (in differences from that of full CI) calculated by the RDM method 
imposing the P,Q,G,T1, or P,Q,G,T1,T2' conditions, and those obtained by CCSD(T) for the 
57 atomic and molecular systems of [45] and the 5 of Table II. 



24 



0.00279 Hartree to full CI for CCSD(T) is for the Oxygen atom at the l D state [45]. 



VI. CONCLUDING REMARKS 

In this chapter, we showed the outline of the reduced-density-matrix method with ap- 
plications to atomic and molecular fermionic systems. Some feature of this method are: 
(i) it is an ab initio method, which is rigorously the same as the Schrodinger equation for 
the ground state; (ii) the number of variables is always four, regardless of the size of the 
system; (iii) from the sparsity of the first- and second- order reduced density matrices the 
existence of a linear scaling method is apparent. The major obstacle for this method is 
the fundamentally difficulty of obtaining the complete iV-representability conditions for the 
2-RDM. However we know fairly good approximated (necessary) condition like P, Q, G, Tl 
and T2' conditions which reproduces comparable ground state energies to CCSD(T), which 
is considered the golden standard method in quantum chemistry. The considered problem 
becomes a semidefinite programming problem, which minimizes a linear functional keeping 
the eigenvalues of matrices non-negative. In this chapter, we presented new results; NFLJ, 
Cn 2 , NH 3 , C2, CH 3 , and using a supercomputer with a highly efficient semidefinite 
programming solver, SDPARA. The semidefinite programming problems for NH 3 , C2, CH 3 , 
and are the largest problems solved so far in the standard formulation. The correlation 
energies using P, Q, G, Tl, T2f for NH~, CH 2 , NH 3 was 100.4%, for CH 3 was 100.3%, and 
for C 2 was 101.2%, respectively. For we used the double-C basis, and attained 132% of 
correlation energy since the open-shell systems are difficult and due to large space, we only 
employ the P, Q and G conditions. 

We would like to close this chapter saying that the RDM method is a promising method for 
quantum chemistry or condensed matter physics. Developing the RDM method is important 
and fundamental to chemistry and physics. 
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